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Abstract. We present radiation hydrodynamics simulations of the collapse of massive pre-stellar 
cores. We treat frequency dependent radiative feedback from stellar evolution and accretion 
luminosity at a numerical resolution down to 1.27 AU. In the 2D approximation of axially 
symmetric simulations, it is possible for the first time to simulate the whole accretion phase of 
several 10 yr for the forming massive star and to perform a comprehensive scan of the parameter 
space. Our simulation series show evidently the necessity to incorporate the dust sublimation 
front to preserve the high shielding property of massive accretion disks. Our disk accretion 
models show a persistent high anisotropy of the corresponding thermal radiation field, yielding to 
the growth of the highest-mass stars ever formed in multi-dimensional radiation hydrodynamics 
simulations. Non-axially symmetric effects are not necessary to sustain accretion. The radiation 
pressure launches a stable bipolar outflow, which grows in angle with time as presumed from 
observations. For an initial mass of the pre-stellar host core of 60, 120, 240, and 480 M© the 
masses of the final stars formed in our simulations add up to 28.2, 56.5, 92.6, and at least 
137.2 M Q respectively. 
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1. Introduction 



The understanding of massive stars suffers from the lack of a generally accepted for- 
mation scenario. If the formation of high-mass stars is treated as a scaled-up version of 
low-mass star formation, a special feature of these high-mass proto-stars is the inter- 
action of the accretion flow with the strong irradiation e mitted by the ne wborn stars 
due to their short Kelvin -Helmholtz contraction timescale ( Shu et al. 1987). Early one- 
dimensional studies (e.g. Larson fc Starrfield 1971 ; Kahnl [l974 : Yorke fc Krugel 1977 ) 
agree on the fact that the growing radiation pressure potentially stops and reverts the 
accretion flow onto a massive star. B ut this radiati ve impact strongly depends on the 
geometry of the stellar environment ( Nakanol 1989t ). The possibility was suggested to 
overcome this radiation pressure barrier via the formation of a long-living massive cir- 
cumstellar disk, which forces the genera tion of a strong aniso t ropic feature of the thermal 
radiation field. Earlier investigations by I Yorke &: Sonnhalter ( 2002 ) tried to identify such 
an anisotropy, but their simulations show an early end of the disk accretion phase shortly 
after its formation due to strong radiation pressure feedback. 
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2. Method 

2.1. The aim 

Aim of this study is to reveal the details of the radiation-dust interaction in the vicinity 
of the most massive star formed during the collapse of a pre-stellar core (gravitationally 
unstable fragment of a molecular cloud). Due to this strong focus onto the core center, 
the physics in the outer core region have to be simplified, e.g. the usage of quiescent 
initial conditions suppresses further fragmentation of the pre-stellar core. 

2.2. The code 

For this purpose, we use our newly developed self-gravity radiation hydrodynamics code. 
The evolution of the gas density, velocity, pressure, a nd total energy dens ity is com- 
puted using the magneto-hydrodynamics code Pluto3 ( Mignone et all [2007;. including 



full tensor viscosity. The derivation and numerical details of the newly developed fre- 
quency dependent hybrid radiation transport method are summarized by iKuiper et al.l 



(2010a|). Our implementation of Poisson 's equation as well as the description of the dust 



and stellar evolution model are given in lKuiper et al.l (|2010bf ) 



The simulations are performed on a time independent grid in spherical coordinates. 
The radially inner and outer boundary of the computational domain are semi-permeable 
walls, i.e. the gas can leave but not enter the domain. The resolution of the non-uniform 
grid is chosen to be (Ar x r A9) min = 1.27 AU x 1.04 AU around the forming massive 
star and decreases logarithmically in the radial outward direction. The accurate size of 
the inner sink cell is determined in a parameter scan presented in Sect. 13.11 

2.3. Initial Conditions 
Our basic initial condition is very similar to the one used bv lYorke fc Sonnhalter ( 2002 ). 



We start from a cold (T = 20 K) pre-stellar core of gas and dust. The initial dust to gas 
mass ratio is chosen to be Mdust/M gas — 1%. The model describes a so-called quiescent 
collapse scenario without turbulent motion (u r = ug = 0). The core is initially in slow 



r 



rigid rotation (Ju^l/i? = Slo = 5 * 10 Hz). The initial density slope drops with 
and the outer radius of the cores is fixed to r max = 0.1 pc. The total mass M core varies 
in the simulations from 60 up to 480 M Q . 



3. Results 

3.1. The dust sublimation front 

In the following, we check the dependency of the stellar accretion rate on the radius r m i n 
of the inner sink cell in four simulations with r min = 1, 5, 10, and 80 AU. We follow the 
long-term evolution of the runs for at least 10 5 yrs. The results are displayed in Fig. [TJ 

The chosen location r mul of the inner boundary of the computational domain influences 
the resulting accretion rate in two distinguishable effects: First, a smaller sink cell leads 
to a shorter free fall epoch, i.e. to an earlier onset of the disk formation phase. This effect 
of the artificial inner cut-off of the gas disk results on one hand in an overestimation of 
the final mass of the central star by approximately 1 Mq or below (left panel of Fig. [I}, 
but on the other hand does not influence the disk accretion epoch (right panel of Fig. [1]). 

The second effect is related to a potential cut-off of the inner dust disk. The region 
in the vicinity of the forming massive star will be heated up to temperatures beyond 
the dust sublimation temperature and a gap is formed between the central star and the 
dust disk. For an inner sink cell radius r m i n larger than the dust sublimation radius r su bi, 
the region of radiative feedback is artificially shifted to higher radii including a strong 
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decrease in density, opacity, and gravity. As a result, the dust disk looses its shielding 
property, the thermal radiation field retains in major parts its isotropic character and 
the radiation pressure therefore stops the emerging disk accretion phase (cp. the case 
of r m i n = 80 AU in Fig. [TJ. This dependency of the radiation pressure on the radius of 
the sink cell explain s also the abrupt end of the accretion phase in the simulations by 
Yorke fc Sonnhalter (2002), who presented simulations of collapsing pre-stellar cores of 



M COIe = 30 Mq, 60 Mq and 120 M© with a radius of the inner sink cell chosen to be 40, 
80, and 160 AU respectively. Our subsequent simulations meet this concern by using an 
adequate central sink cell radius of r m i n = 10 AU. 

3.2. The radiation pressure barrier 

The spherically symmetric accretion flow simulations yield a maximum stellar mass of 
less than 40 M n independent of the initial core mass M corc > 60 M Q due to radiative 
feedback ([Kuiper et al. I l2010bh . We attack this radiation pressure barrier in axially and 



midplane symmetric circumstellar disk geometry now. Resulting accretion histories as a 
function of the actual stellar mass for different initial core masses of M cole = 60 Mq, 
120 M Q , 240 M Q , and 480 M Q are displayed in Fig. El In face of the additional centrifugal 
forces, the disk accretion easily breaks through the upper mass limit of the final star of 
Af* < 40 M© found in the spherically symmetric accretion models! 
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Figure 2: Accretion rate M* as a function of the actual stellar mass M* for four different 
initial core masses. Also the periods of accretion are mentioned for each run. The vertical 
line marks the upper mass limit found in the spherically symmetric accretion models. 
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No upper mass limit of the final star is detected so far, but the star formation efficiency 
declines for higher mass cores. The depletion of the envelope by radiative forces decreases 
the large-scale accretion onto the midplane. The disk looses its shielding property and 
the radiation pressure starts to accelerate the remnant material in the outward direction. 



4. Summary 

We performed high-resolution radiation hydrodynamics simulations of monolithic pre- 
stellar core collapses including frequency dependent radiative feedback. The dust subli- 
mation front of the forming star could be resolved down to 1.27 AU. The whole accretion 
phase of several 10 5 yrs was computed. The frequency dependent ray-tracing of our newly 
developed radiation module denotes the most realistic radiation transport method used 
in multi-dimensional hydrodynamic simulations of massive star formation by now. 

The broad parameter studies, especially regarding the size of the sink cell and the initial 
core mass, reveal new insights of the radiative feedback onto the accretion flow during the 
formation of a massive star: In the case of disk accretion, the thermal radiation field gen- 
erates a strong anisotropic feature. We found that it is strict necessary to include the dust 
sublimation front in the computational domain to reveal the persistent anisotropy during 
the long-term evol ution of the accretion disk. The short accretion phases of the disks in 
the simulations bv lYorke fc Sonnhalterl f|2002h are a result of the fact that they did noi; 
include the dust sublimation front in their simulations, as clearly shown in our result of 
the parameter scan of the size of the central sink cell (see Sect. l3.lj ). Additional feeding of 
the disk by unstable outflo w regions due to t he so- called "3D radiative Rayleigh- Taylor 
instability" , as proposed in iKrumholz et~aT1 fj2009h . would enhance this anisotropy but 
is not necessary. In fact, preliminary analyses of our ongoing three-dimensional collapse 
simulations identify evolving gravitational torques in the massive circumstellar disk as 
well as the launching of non-axially symmetric outflows. The accretion rate, driven by the 
angular momentum transport of the gravitational torques, is more episodically compared 
to the axially-symmetric runs, but results in a similar mean accretion rate. The radiation 
pressure driven ouflow remains stable. 

Finally, the central stars in our simulations of the disk accretion scenario grow far be- 
yond the upper mass limit found in the case of spherically symmetric accretion flows. For 
an initial mass of the pre-stellar host core of 60, 120, 240, and 480 M Q the masses of the 
final stars formed add up to 28.2, 56.5, 92.6, and at least 137.2 M respectively. Indeed, 
the final massive stars are the most massive stars ever formed in a multi-dimensional 
radiation hydrodynamics simulation so far. 
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